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METHOD OF DIGITAL IMAGE ENHANCEMENT AND SHAB PENING^ 

Priority Claim: This application claims the benefit of the filing date of United States 
Provisional Patent Application Serial Number 60/100,136, filed September 14, 1998 for 
"METHOD OF DIGITAL IMAGE ENHANCEMENT AND SHARPENING BASED ON 
MINIMUM GRADIENT SUPPORT CONSTRAINT." 

Technical Field: The present invention relates in general to image processing, and, in 
particular, to image enhancement and sharpening. 

The method, for example, can be applied to optical image processing, for image restora- 
tion and sharpening in biomedical, geophysical, astronomical, high definition television, re- 
mote sensing, and other industrial applications. 

Background Art: Comprehensive coverage of the prior art may be found, for example, 
in W. K. Pratt, "Digital Image Processing," 2nd Edition, John Wiley and Sons, NY (1988); 
H. Stark, "Image Recovery; Theory and Application," Academic Press, Inc., Harcout Brace 
Jovanovich Publishers, New York (1987); R. C. Gonzalez and P. Wintz, "Digital Image Pro- 
cessing," 2nd Edition, Addison- Wesley Publishing Company, Inc., Advanced Book Program, 
Reading, Mass. (1987); and R. L. Lagendijk and J. Biemond, "Iterative Identification and 
Restoration of Images," Kluwer International Series in Engineering and Computer Science, 
Kluwer Academic Publishers, Boston, Mass., (1991). 

There have been several attempts to develop a method of image processing and restoration 
based on the solution of the linear ill-posed inverse problem: 

d = Bm, (1) 

where d is the blurred (or degraded) image, m is original (or ideal) image, and B is the 
blurring linear operator of the imaging system. Note that the original image, as well as the 
blurred image, can be defined in a plane (2-D image: m = m(x,y) , d = d(x,y)) or in a 
volume (3-D image: m = m(x,y,z) , d = d (x, y, z)) . 

A wide variety of electron-optical devices obey equation (1) with different blurring op- 
erators as noted by C. B. Johnson et al., in "High-Resolution MicroChannel Plate Image 
Tube Development," Electron Image Tubes and Image Intensifies II, Proceedings of the 
Society of Photo-Optical Instrumentation Engineers, Vol. 1449. I. P. Csorba, Ed. (1991), 
pp. 2-12. These devices are used in various biomedical imaging apparatus, including image 
intensifier-video camera (II-TV) fluoroscopic systems (see S. Rudin et al., "Improving Fluo- 
roscopic Image Quality with Continuously Variable Zoom Magnification," Medical Physics, 
Vol. 19 (1991), PP- 972-977); radiographic film digitizers (see F. F. Yin et al., "Measurement 
of the Presampling Transfer Function of Film Digitizers Using a Curve Fitting Technique," 
Medical Physics, Vol. 17 (1990), pp. 962-966); radiographic selenium imaging plates (see P. 
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J. Papin and H. K. Huang, "A Prototype Amorphous Selenium Imaging Plate System for 
Digital Radiography," Medical Physics, Vol. 14 (1987), pp. 322-329); computed radiography 
systems (see S. Sanada et al.. "Comparison of Imaging Properties of a Computed Radio- 
graphy System and Screen-Film Systems," Medical Physics, Vol. 18 (1991), pp. 414-420; 
H. Fujita et al., "A Simple Method for Determining the Modulation Transfer Function in 
Digital Radiography," IEEE Transactions on Medical Imaging, Vol. 11 (1992), pp. 34-39); 
digital medical tomography systems (see M. Takahashi et al., "Digital 'IV Tomography: 
Description and Physical Assessment," Medical Physics, Vol. 17 (1990), pp. 681-685). 

Geophysical, airborne, remote sensing, and astronomical blurred images also can be de- 
scribed by equation (1) (see M. Bath, "Modern Spectral Analysis with Geophysical Applica- 
tions," Society of Exploration Geophysicists, (1995), 530 pp.; C. A. Legg, "Remote sensing 
and geographic information systems," John Wiley & Sons, Chichester, (1994), 157 pp.). 

Most prior efforts to solve the problem: (1) are based on the methods of linear in- 
verse problem solutions. Inverse problem ( 1) is usually ill-posed, i.e., the solution can be 
non-unique and unstable. The conventional way of solving ill-posed inverse problems, ac- 
cording to regularization theory (A. N., Tikhonov, and V. Y., Arsenin, "Solution of ill-posed 
problems," V.H. Winston and Sons., (1977); M.S., Zhdanov, "Tutorial: regularization in in- 
version theory," Colorado School of Mines (1993)), is based on minimization of the Tikhonov 
parametric functional: 

P a {m) = <f>{m)+as{m), (2) 

where 4> ( m ) * s a niisfit functional determined as a norm of the difference between observed 
and predicted (theoretical) degraded images: 

if (m) = \\Bm - df = (Bm - d, Bm - d) . (3) 

Functional s (m) is a stabilizing functional (a stabilizer). 

It is known in the prior art that there are several common choices for stabilizers. One is 
based on the least squares criterion, or, in other words, on L 2 norm for functions describing 
the image: 

■?L 2 (m) = |fm.|| 2 = (m,m) = j m 2 dv = min, (4) 

where V is the domain (in 2-D space or in 3-D space) of image definition, and (..., ...) denotes 
the inner product operation. 

The conventional argument in support of this norm comes from statistics and is based on 
an assumption that the least square image is the best over the entire ensemble of all possible 
images. — 

Another stabilizer uses minimum norm of difference between the selected image and some 
a priori image m apr : 

SL 2 «pr<m) = \\ m ~ m-aprf = min . (5) 
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This criterion, as applied to the gradient of image parameters Vm, brings us to a maximum 
smoothness stabilizing functional: 

SmaxsTn (m) = |j Vm|| 2 = (Vm, Vm) = min . (6) 

Such a functional is usually used in inversion schemes. This stabilizer produces smooth 
images. However, in many practical situations the resulting images don't describe properly 
the original (ideal) image. Inversion schemes incorporating the aformentioned functional 
also can result in spurious oscillations when m is discontinuous. 

It should be noted that in the context of this disclosure, "image parameters" is intended 
to describe the physical properties of an examined media. Such parameters include without 
limitation: color, intensity, density, field strength, and temperature. Discrete samples of 
such parameters are assembled, eg. in a matrix [dj, to describe an image of the media. 

In the paper by (L.I., Rudin, Osher, S. and Fatemi, E., "Nonlinear total variation based 
noise removal algorithms," Physica D, 60, (1992), pp. 259 - 268) a total variation (TV) based 
approach to reconstruction of noisy, blurred images has been introduced. This approach uses 
a total variation stabilizing functional, which is essentially L\ norm of the gradient: 

stv (m) = \\Vm\\ Ll = j y |Vm| dv. (7) 

This criterion requires that an image's parameter distribution in some domain V be of 
bounded variation (for definition and background see (E., Giusti, "Minimal surfaces and 
functions of bounded variations," Birkhauser (1984)). However, this functional is not differ- 
entiable at zero. To avoid this difficulty, Acar and Vogel introduced a modified TV stabilizing 
functional (R., Acar, and C.R., Vogel, "Analysis of total variation penalty methods," Inverse 
Problems, 10, (1994), pp. 1217-1229): 

sprv (m) = J \f\Vm\ 2 + e 2 dv. (8) 

The advantage of this functional is that it doesn't require that the function m is continuous, 
but just that it is piecewise smooth (C.R., Vogel, and M.E., Oman, "A fast, robust algo- 
rithm for total variation based reconstruction of noisy, blurred images," IEEE Transactions 
on Image Processing, (1997)). The TV norm doesn't penalize discontinuity in the image pa- 
rameters, so we can remove oscillations while sharp conductivity contrasts will be preserved. 
At the same time it imposes a limit on the total variation of m and on the combined arc 
length of the curves along which m is discontinuous. 

TV functionals stv (m) and sprv (tti), however, tend to decrease bounds of the image 
parameters variation, as can be seen from (7) and (8), and in this way still try to "smooth" the 
real image, but this "smoothness" is much weaker than in the case of traditional stabilizers 
(6) and (5). 

In yet another development A. S. Carasso introduced a procedure for digital image 
restoration based on Tikhonov regularization theory and description of the blurred image d 
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as the end result of a diffusion process appiied to the desired ideal image m (A.S., Carasso, 
"Overcoming Holder Continuity in 111- Posed Continuation Problems," SIAM Journal on 
Numerical Analysis, Volume 31, No. 6, December 1994, pp. 1535-1557; C.E., Carasso, 
United States Patent, 5,627,918, 5/1997). In the framework of Carasso's method the image 
restoration procedure is equivalent to solving diffusion equation backwards in time using the 
given degraded image d as data at the final moment of diffusion. The desired ideal image is 
the corresponding solution of the diffusion equation at the initial time moment. To generate 
a stable solution, Carasso suggests to use a constraint 

||m - B $ m\\ = min, 

where B s is the fractional power of the operator B. 

The main limitation of Carasso's approach is that it can be applied only to a specific class 
of integral blurring operators, which can be described by the diffusion problem solution. 

In yet another development (H. Kwan, and J.T., Liang, Digital image sharpening method 
using SVD block transform. United States Patent, 4,945,502, 7/1990), the authors suggest 
to use the Singular Value Decomposition technique and a priori known information about 
the statistical properties of the noise and image texture to suppress the noise in the degraded 
images. 

The foregoing attempts met with varying degrees of success in developing the methods 
of digital image restoration. However, from the point of view of practical applications, these 
methods are not entirely satisfactory. Therefore, it is highly desirable to develop a method 
which can be applied to a wide variety of image restoration problems. 

DISCLOSURE OF THE INVENTION 

In the method of the present invention, a new approach to digital image enhancement and 
sharpening is developed based on a new type of constraint on the solution which generates 
focused and sharp images. 

The method of the present invention may be used in industrial applications, including 
biomedical imaging for human body study, for example, in breast cancer diagnosis, (see S. 
Webbet et ah, "Constrained Deconvolution of SPECT Liver Tomograms by Direct Digital 
Image Restoration," Medical Physics, Vol. 12 (1985), pp. 53-58; U. Raft et al, "Im- 
provement of Lesion Detection in Scintigraphic Images by SVD Techniques for Resolution 
Recovery," IEEE Transactions on Medical Imaging, Vol. MI-5 (1986), pp. 35-44; B.C. 
Penney et al., "Constrained Least Squares Restoration of Nuclear Medicine images; Select- 
ing the Coarseness Function," Medical Physics, Vol. 14 (1987), pp. 849-859; B.C. Penney 
et al., "Relative Importance of the Error Sources in Wiener Restoration of Scintigrams," 
IEEE Transactions on Medical Imaging, Vol. 9 (1990), pp. 60-70; K. S. Pentlow et al., 
"Quantitative Imaging Using Positron Emission Tomography with Applications to Radioim- 
munodiagnosis and Radioimmunot! herapy," Medical Physics, Vol. 18 (1991), pp. 357-366); 
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magnetic resonance imaging (see S. M. Mohapatva et al., "Transfer Function Measurement 
and Analysis for a Magnetic Resonance Imager," Medical Physics, Vol. 18 (1991), pp. 1141- 
1144); and computed tomography scanners (see E. L. Nickoloff and R. Riley, "A Simplified 
Approach for Modulation Transfer Function Determinations in Computed Tomography," 
Medical Physics, Vol. 12 (1985), pp. 437-442). 

The method of the present invention may also find application in night vision systems; 
imaging through the atmosphere; in observational astronomy (see N. S. Kopeika, "Imaging 
Through the Atmosphere for Airborne Reconnaissance," Optical Engineering, Vol. 26 (1987), 
pp. 1146-1154; J. D. Gonglewski and D. G. Voelz, "Laboratory and Field Results in Low 
Light Post-detection Turbulence Compensation Using Self Referenced Speckle Holography," 
Digital Image Synthesis and Inverse Optics, Proceedings of the Society of Photo-Optical 
Instrumentation Engineers, Vol. 1351, A. E Gmitro, P.S. Idell, and I. J. LaHaie, Eds. 
(1990), pp. 798-806). 

The instant method may also find excellent application in geophysical exploration for 
generating sharp resolved seismic, electromagnetic and gravity images of underground geo- 
logical structures, for undersea imaging, in airborne geological and geophysical sensing and 
remote sensing for imaging from satellites (see M. Bath, "Modern Spectral Analysis with 
Geophysical Applications," Society of Exploration Geophysicists, (1995), 530 pp.; C. A. 
Legg, "Remote sensing and geographic information systems," John Wiley & Sons, Chich- 
ester, (1994), 157 pp.). 

Other important image sharpening applications of the present method include high defi- 
nition television (HDTV) (see N. V. Rao, "Development of a High-Resolution Camera Tube 
for 2000- Line TV System," Electron Tubes and Image Intensifies, Proceedings of the So- 
ciety for Photo-Optical Instrumentation Engineers, Vol. 1243, 1. P. Csorba, Ed. (1990), pp. 
81-86). 

The method of the present invention introduces a new approach to reconstruction of 
blurred images which is based on using a specially selected constraint or stabilizing functional 
(stabilizer). In particular, a new form of a constraint (stabilizer) which minimizes the area 
where strong image parameters variations and discontinuity occurs is developed. This new 
constraint may be characterized as a minimum gradient support (MGS) constraint. A MGS 
constraint in combination with a penalization function helps to generate a stable solution 
of the inverse linear problem ( 1) for arbitrary blurring imaging operator. The procedure of 
the instant method finds practical application to digital image enhancement in general. 

The present invention provides a new method for digital image enhancement and sharp- 
ening in a system which can be described by the linear operator equation in a form d = 
Bm, where d is degraded (blurred) image, m is original (ideal) image, and B is the blur- 
ring operator of the imaging system. In this method, the norm of difference between the 
observed degraded image d and numerically predicted degraded image d pr corresponding 
to the sharp ideal image m s is preassigned to be equal to a given tolerance value e : 
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— dpr\\ = \\d— Bm s \\ = e. where denotes L 2 norm. The restored (sharpened) im- 
age is selected from the class of possible ideal images by minimization of the quantity: 
- — H^lHL — =■ — minimum, under the constraint \\d~ Bm s \\ — £, where Vm s describes the 
spatial gradient of the image parameters, and e C 1 is a substantially small constant. 
The method can be implemented numerically using a general purpose computer, or it may 
be performed using dedicated hardware specially designed to solve constraint minimization 
problems. 



BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 is a flow chart illustrating the image enhancement and sharpening method of this 
invention. 

FIG. 2 is an ideal seismic image of a geological cross-section. 

FIG. 3 is a typical blurred seismic image of the same cross-section shown in Figure 2 
which is generated as the result of conventional seismic exploration technique. 

FIG. 4 is a reconstructed image generated by the method of digital image enhancement 
and sharpening of the present invention. 

BEST MODE FOR CARRYING OUT THE INVENTION 

The principles of image enhancement and sharpening can be described as follows. 

We introduce a new constraint (stabilizer) which minimizes the area where strong image 
parameters variations and discontinuity occur. "We call this new constraint (stabilizer) a 
minimum gradient support (MGS) constraint. We demonstrate that a MGS constraint in 
combination with the penalization function helps to generate a stable solution of the linear 
inverse problem ( 1) describing the image processing. We call this approach sharpening 
blurred images. 

For the sake of simplicity we will discuss first a minimum support (MS) functional, which 
provides the image with the minimum area of the anomalous image parameters distribution. 
Consider the following integral of the image parameters distribution in domain V: 

We introduce the support of m (denoted sptm) as the combined closed subdomains of V 
where m ^ 0. We call sptm an image parameter support. Then expression (9) can be 
modified: 2 ^ 

J ' {rn) = C- 1 1 ~ dv = sp ' m ~ e2 4t„ ^T?*'- (10) 

i,From the last expression (10) we can see that: 



J e (m) — > sptm, if e — > 0. 



(ID 
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Thus, integral J e (m) can be treated as a functional, proportional (for a small e) to the image 
parameter support. We can use this integral to introduce a minimum support stabilizing 
functional Spus (m) as: 

(m - m apT ) 2 



s (m) = J e (m~ m apr ) = f {m m ^> dv (12) 
Jv (m - m apr y + e 2 v J 



To justify this choice one can prove that s ms (m) can actually be considered as a stabilizer 
according to regularization theory (A. N., Tikhonov, and V. Y., Arsenin, "Solution of ill- 
posed problems," V.H. Winston and Sons). 

This functional has an important property: it minimizes the total area with nonzero 
departure of the image parameters from the given a priori image m apr . Thus, a dispersed 
and smoothed distribution of the degraded image parameters with all values different from 
the a priori image m apr results in a big penalty function, while well-focused distribution 
with a small departure from m apr will have a small penalty function. A similar approach 
was originally used in the paper (B.J., Last, and K., Kubik, "Compact gravity inversion," 
Geophysics, 48, (1983), pp. 713-721) for compact gravity inversion. 

We can use this property of a minimum support functional to increase the resolution of 
imaging system. To do so we modify s eMS (m) and introduce a minimum gradient support 
junctional as: 

SeMGS (m) = J e [Vm] = f Vffl _ Vm dv. (13) 
JvVm-Vm + e 2 K J 

We denote by sptVm the combined closed subdomains of V where Vm ^ 0. We call 
sptVm an image parameter gradient support (or, simply, gradient support). Then expression 
(13 ) can be modified: 



^mgs (m) = sptVm -e 2 f — ~ dv. 

/sptVm Vm ■ Vm 4- e 2 



/sptVm Vm ■ Vm 4- e 2 

iFrom the last expression we can see that: 



(14) 



SeMGS (m) -> sptVm, if e 0. (15) 

Thus, we can see that s eMGS (m) can be treated as a functional, proportional (for a small e) 
to the gradient support. One can demonstrate that the minimum gradient support functional 
satisfies Tikhonov criterion for a stabilizer. 

Thus we arrive at a new method of the present invention for blurred image restoration. In 
this method the norm of difference between the observed degraded image d and numerically 
predicted degraded image d^ corresponding to the sharp ideal image m s is preassigned to 
be equal to a given tolerance value £ : 



Wd-d^W^Wd-BmsW^e, 



(16) 
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where ||...|| denotes L 2 norm. The restored (sharpened) image is selected from the class of 
possible ideal images by minimization of the quantity: 

l!Vm s || 2 

tt= = 7, r = minimun, ( 17\ 

||Vm s • Vm s || + e 2 1 ; 

under the constraint 

\\d- Bm s \\ = e, (18) 

where Vm s describes the spatial gradient of the image parameters, and e < 1 is a small 
constant. 

It can be demonstrated that the solution of the minimization problem (17) under the 
constraint (18) is reduced to the minimization of the unconstrained Tikhonov parametric 
functional P a (equation (2)), which can be written as follows: 

P a (m) = (Bm ~d,Bm-d)+ as eM Gs{rn), (19) 

where(..., ...) denotes the inner product operation determined in (4), and a is regulaxization 
parameter. 

The regularization parameter a describes the trade-off between the best fitting and rea- 
sonable stabilization. In a case when a is selected to be too small, the minimization of 
: the parametric functional P a (m) is equivalent to the minimization of the misfit functional 
= <f>{m) = {Bm - d, Bm - d) =minimum, therefore we have no regularization, which can result 
in an unstable incorrect solution. When a is too large, the minimization of the parametric 
functional P a (m) is equivalent to the minimization of the stabilizing functional s eM Gs{m), 
which will force the solution to be closer to the a priori model. Ultimately, we would expect 
the final model to be exactly like the a priori model, while the observed data are totally 
ignored in the inversion. Thus, the critical question in the regularized solution of the inverse 
problem is the selection of the optimal regularization parameter a. The basic principles use! 
d for determining the regularization parameter a are discussed in (A. N., Tikhonov, and V. 
Y., Arsenin, "Solution of ill-posed problems," V.H. Winston and Sons (1977)). We suggest 
to use a simple numerical method to determine the parameter a . Consider for example the 
progression of numbers 

a k = a 0 g k ; k = 0, 1, 2, n; q > 0. (20) 
For any number a k we can find an element m ak , minimizing P ak (m), and calculate the misfit 
(f}(m ak ). The optimal value of the parameter a is the number a fc0 , for which we have 

0(«O = £, (21) 
where e a given tolerance value. The equality (21) is called the misfit condition. 

The important element of this invention is that we introduce a new parametrization m g 
of the reconstructed image m, equal to the gradient of original image parametrization: 



m g = Vm. 



(22) 
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We cail m g image gradient parameters. Using the new parametrization m g we can rewrite 
the original equation of the imaging system (1 ) in a form: 

d = B g m 9 , (23) 

where new blurring operator B g is composed of the original blurring operator B and inverse 
gradient (integral) operator V" 1 : 

B g = BV~ l . (24) 

Taking into account notation (22) the MGS functional introduced above can be written 
as the squared L 2 norm of some functions of the model parameters: 



SeMGsim) = (f (m g ) , / (m g )) , (25) 

(iKir 



/(m * )= /„ IT' 1 (26) 



and e is a substantially small number. 

In the method of this invention we introduce a variable weighting function: 

1 

(iKIf 

Then the MGS functional in general case can be written as 



W *^= ,„ , i2 . „m/2 - (27) 



SeMGsim) = (f {m g ) , f (m f )) « (w e (m g ) m g , w e (m g ) m g ) = 
= (m g , m g ) Wt = J w\ (m g ) m 2 g dv 

where (rn g ,Tn g ) Wc denotes the weighted square norm of m g . 
The corresponding parametric functional has the form: 

P a (m) = P«(m g ) = (B g m 9 - d, B g m g - d) + a (m s , m 9 ) Wc . (28) 

Therefore, the problem of minimization of the parametric functional introduced by equa- 
tion (19) can be treated in a similar way to the minimization of the conventional Tikhonov 
functional with the L 2 norm stabilizer s L30pr (m) (equation (5)). The main difference is that 
now we introduce some a priori variable weighting functions w e (m g ) for model parameters. 

Another important element of the present invention is that the computational procedure 
to minimize the parametric functional (28) is based on the re-weighted conjugate gradient 
method in the space of weighted image parameters. 
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Consider a discrete inverse problem equivalent to the linear inverse problem ( 1) in the 
case of the discrete image parameters. Suppose that original and degraded images can be 
represented in discrete form by the vectors m and d. 

In this case, the equation (1) can be rewritten in a matrix notation: 

d = Bm, (29) 

where B is the matrix of the blurring operator. Introducing new parametrization of the 
image m g we can write the matrix equation (29) in a form 

d = B g m g , 

where new blurring operator matrix B g is composed of the original blurring operator matrix 
B and the matrix of the inverse gradient (integral) operator V -1 : 

B 9 = BV~ l . 

The matrix V -1 may be easily precomputed as the inverse matrix to the matrix of the 
finite-difference gradient operation applied to the discrete image parameters. 
Parametric functional (28) can be expressed now using matrix notations: 

P g a (m) = (B g (m 9 ) - dy(B 9 (m a ) - d)+ 

a(W e m g yW e m g , (30) 
where asterisk "*" denotes transposed complex conjugate matrix, and W e — W e (m g ) is the 
matrix of the weighting function w e (m g ) introduced above in equation ( 27). We call matrix 
W t a sharpening weighting matrix. 

The problem of minimization of the parametric functional (30) can be reformulated using 
the space of weighted parameters: 

= W e m g (31) 

We can consider the modified blurring operator B™ which relates the new weighted param- 
eters nig" to the data 

d = B™ (m%) (32) 
In order to keep the same data we should assume 

B? = B 9 W? (33) 

Using these notations we can rewrite the parametric functional as follows 

/?«) = {B g {W?m») - d)* (B g (W e ~^) - d) + 

a{rh%¥fh%. (34) 
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In other words, we keep the same misfit, as in equation (30) because 

<Pu, K) = (B g (W- l m™) - dy(B g (W- l m%) - d) = 

. ~ . ( 35 ) 

{B s {fh s )-dy{B g {m 9 )-d) 

and the same stabilizer, as in equation (30), equal to the least square norm of m g with the 
weights 

Sw = (mj)*mj = rhr g Wlm g . (36) 
Therefore, the minimization problem (34) is equivalent to the minimization problem (30). 
However, the iterations generated in the framework of gradient type methods would be 
different for (34) in comparison with the iterations introduced above for (30). Moreover, 
numerical examples show that, as a rule, iterative process converges faster for (34) than for 
(30). 

To construct the iterative process in the space of weighted parameters we perturb now 
the parametric functional with respect to the weighted parameters fh w : 

6^P« (m-) = 2 (Sn%y \b? (B 9 (m g ) - d) + «nj] (37) 

We can apply first the traditional gradient type method. In this case the perturbation 
of the model parameters should be selected to be equal to: 

5m™ = -kl{m w g ) ( 38 ) 

where l(m g ) is the direction of steepest ascent, 

l(m%) = W-'B* g (B g (W e -*rn$) - d) + am™ (39) 

and k is the length of the step. 

In the method of this invention we propose to use the conjugate gradient method to find 
the minimum of the parametric functional. It is based on the successive line search in the 
conjugate gradient direction l(m n ): 

™£n+i = ™£ n + 5fff° = ro™ n - M(m£ n ). (40) 

The idea of the line search is the following. We present P° (m£ n - M(m£j) as a function 
of one variable k, and, evaluating it three times along direction Z(m£J, approximately fit it 
by a parabola, find its minimum and the value of k n , corresponding to this minimum. 

The conjugate gradient directions /(m£„) are selected as follows. On the first step we 
use the gradient direction: 

^o) = *K„) = W^§; (B 9 (W- 0 l ml Q ) -d)+ ami, (41) 
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On the next step, the conjugate gradient direction is the linear combination of the gra- 
dient on this step and the direction l{m%) on the previous step: 

On the nth step 

= T(n%* + i) + /W(m£ n ). 

where 

KKn) = W-'B; (B 9 (W- n l ml n ) -d) + aml n , 

and 

™g,n = W en m gtn . 

The coefficients p n+x are denned from the condition that the directions f(m£ n+1 ) and i(fn^ B ) 
are conjugate: 

Note that the iterative process (40) converges to the minimum of the functional (30). This 
algorithm always converges because on every iteration we apply the parabolic line search. 
We call this algorithm conjugate gradient re-weighted optimization because the sharpening 
weighting matrix W en is updated on every iteration. One can find the formal proof of the 
convergence of this type of optimization technique in (U., Eckhart, "Weber's problem and 
Weiszfeld's algorithm in general spaces," Math. Programming, v. 18, ( 1980) pp. 186-196). 
Note that in the case of linear forward operator B g parametric functional has only one 
local minimum, so the minimization of P° is unique (A. N., Tikhonov, and V. Y., Arsenin, 
"Solution of ill-posed problems," V.H. Winston and Sons., (1977)). 

The advantage of the conjugate gradient re- weighted optimization algorithm is that we do 
not have to know the gradient of f{m) for every iteration, but only its value for corresponding 
model parameters, which is easy to calculate. 

The final step of the described algorithm involves recalculating the real image parameters 
from the weighted solution according to formula: 

Another important new element of this invention is introducing the penalization function 
in the image restoration process. 

Let us assume that the ideal image can be described as a combination of the background 
image, mb g (r) , and the anomalous part of the image, m a (r) , where r is the coordinate 
of some point in domain V. In this situation it is required that there should be only two 
values of the parameter m (r) in the focused image equal to the background value or to the 
anomalous value. However, the geometrical distribution of these values is unknown. We 
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can force the inversion to produce an image which not only generates the observed degraded 
image but which is also described by these known values, thus painting the geometry of the 
ideal object. We call this approach "penalization". There is a simple and straightforward 
way of combining penalization and the MGS method. Numerical tests show that MGS 
generates a stable image, but it tends to produce the smallest possible anomalous domain. 
It also makes the image look u! nrealistically sharp. At the same time, the image parameter 
values m (r) outside this local domain tend to be equal to the background value m bg (r) , 
which nicely reproduces the background image. We can now impose the upper bound for 
the positive anomalous parameter values mf (r) , and during the sharpening process simply 
chop-off, or discard, all the values above this bound. This algorithm can be described by 
the following formula: 

m (r) - m bg (r) = < 6 (r) , if [m (r) - m bg (r)] > m« b (r) ; 
m (r) - m bg (r) = 0, if [to (r) - m bg (r)] ^ 0. 

Thus, according to the last formula the image parameters values m (r) are always distributed 
within the interval: 

m bg (r) <m{v) < rn bg (r) + (r) . 
A similar rule is applied for the case of negative anomalous parameter values. 

Hence, an important feature of the present invention is the ability to enhance and sharpen 
degraded images. It is realized numerically as the following seven steps. 

Step 1. Constructing the initial restoration of the ideal image gradient parameters m gA 
according to the formula 

m sA = -MO = -k 0 W^B* g (B g (W^ fi ) -d)~ akomlo (43) 

Note that we assume that starting zero iteration is equal to zero fffg fi = 0, and W e0 is equal 
to the identity matrix W e0 = I, so the starting approximation of the ideal image gradient 
parameters is generated by applying the transposed complex conjugated blurring operator 
and the inverse gradient operator to the blurred image. Mathematically it can be described 
by matrix multiplication: 

rh 9tl = k 0 B*d= koV^B'd, (44) 
where k 0 is the length of the initial step. 

Step 2. Computing the inverse sharpening filter (matrix W^ 1 = W~ l e (m 9il ) on the 1-st 
iteration as the matrix of the function w e (m ?)1 ) introduced above in equation (27): 

™ e K,i) = (llm^H 2 + e 2 ) 1/2 « j|m Pil j| . (45) 

Step 3. Constructing the partially sharpened image (second iteration to the weighted 
restored image) according to the formula: 
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where 



am 



'9,0' 



The length of the step ki is determined by the line search to the minimum of the functional 
Pg _ k K™%,i)) • Parameter is calculated by the formula: 



||'(Ko)|| ' 

Step 4. Computing the inverse sharpening filter (matrix W en ) on the n-th iteration 
(n = 2, 3, ...) as the matrix of the weighting function w e {m gM -i) introduced above in equation 
(27): 

w; 1 K,^) = (\\m 9 , n -if + e 2 ) 1/2 « Wm^W . (46) 

Step 5. Constructing the partially sharpened weighted image {{n + l)-th iteration to 
the restored image) according to the formula 

™s,n+i = + &™%, n = fn£ n - fcj(m£ n ), (47) 
where weighted conjugate gradient directions are determined by the expressions: 

K^ln) = Krn w 9 ,n) + A.* (m£ n - J, (48) 
K^ln) = W- l B* s (B, (W-^ln) -d)+ aml n . (49) 
The length of the step k n is determined by the line search to the minimum of the functional 
P a (mjj^ - kj(™%, n )) • Parameter 0 n is calculated by the formula: 

0 ii r K-)ir 

Step 6. Recomputing the real parameters of the ideal image from the weighted param- 
eters at the (n + l)-th iteration by inverse filtering the weighted image: 

™ n+ i=V-iW~l +1) ^ n+1 (50) 

Step 7. Imposing the upper bound m% b (r) for the positive anomalous parameter values, 
and the lower bound -m l £ (r) for the negative anomalous parameter values and chopping-off, 
or discarding, all the values outside these bounds during entire image restoration process to 
keep the image parameters values m n (r) always distributed within the interval: 

m bg (r) - (r) X m n (r) ^ m hg (r) + mf (r) . (51) 
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In the drawing, Figure 1 presents the flow chart illustrating the image enhancement and 
sharpening method 10 of the present invention. The method can be performed upon the 
degraded image expressed as a digitized M x M array d = d (arj, yj){fov a 2-D image) or a 
digitized MxM xM array d = d {x it yj, z k ) (for a 3-D image). The linear blurring operator 
B is assumed to be given as a square matrix of the size M 2 x M 2 (for a 2-D image) or 
M 3 x M 3 (for a 3-D image). 

Referring again to Figure 1, the first step is the initial restoration of the image (block 
11) by applying a transposed matrix of blurring operator and a inverse matrix of gradient 
operator (integration) to the digital blurred (degraded) image according to formula (43). 

An inverse sharpening filter, according to the present invention, is then numerically con- 
structed, as represented in block 12 (Figure 1), by use of a function w e (m 9)1 ) accomplishing 
a minimum gradient support constraint according to formula (45). 

In partially sharpened image block 13 of image enhancement and sharpening method 10, 
the subsequent iteration of the weighted restored image is constructed according to formulae 
(47), (48 ), and (49). 

In inverse filtering block 14, the real parameters of the ideal image are recalculated from 
the weighted parameters by inverse filtering the weighted restored image. 

In penalization block 15, the image parameters values are forced to be distributed within 
the given interval to produce a sharp, enhanced original image. 

The construction of the sharpening filter is then undone in block 16. This is the same 
procedure as used in block 12 of image enhancement and sharpening method 10, but applied 
now to the image constructed in the previous block 15 acording to formula (46). 

Execution of the image enhancement and sharpening method 10 may then return in a 
recursive loop through constructing sharpening filter block 12, to partially sharpened image 
block 13, to inverse filtering block 14, and finally to penalization block 15, until the norm 
of a difference between the observed degraded image d and numerically predicted degraded 
image dp r corresponding to the sharp ideal image generated in block 15 reaches a given 
tolerance value e : JJrf - (^.Jj = jjd - Bm^jjj = e. 

An important feature of the present invention is that the image deblurring method 10 
can be applied to restoration of blurred images of arbitrary origin. 

As an example of an application of this method, we present here a result of an image 
enhancement and sharpening experiment conducted for a seismic image. Figure 2 shows the 
ideal seismic image of a geological cross-section. Figure 3 presents the typical blurred seismic 
image of the same cross-section which is generated as the result of conventional seismic 
exploration technique. This image has been processed by the digital image enhancement 
and sharpening method 10 of the present invention. The reconstructed image is shown in 
Figure 4. One can see that this image is practically identical to the original ideal image 
presented in Figure 2. 

Hence, an important feature of the present invention is the ability to produce a real- 
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istic original image from a degraded one. It is believed that the image enhancement and 
sharpening method 10 of the present invention can be used in optical image processing, for 
image restoration and sharpening in biomedical, geophysical, astronomical, high definition 
television, remote sensing, and other industrial applications. 

The method of the present invention is based on three new ideas that make it useful 
in image restoration technique. The first idea is to use a new physical constraint in image 
restoration procedure based on minimization of the area where strong image parameter varia- 
tions and discontinuities occur. This constraint provides substantial qualitative improvement 
in image reconstruction procedures. The second idea is to implement this constraint in a 
form of weights imposed on image parameters, which makes it practical to realize this con- 
straint in a numerical algorithm. The third idea is penalization of the image parameters to 
keep them within given reasonable limits of variations. Together these three new ideas result 
in a new method 10 of digital image enhancement and sharpening. 

The method 10 of the present invention can be implemented numerically using a general 
purpose computer programmed to perform the procedures of the blocks 11-16. It can 
be understood that the method 10 of the present invention may also be performed using 
dedicated hardware specially designed to solve the constraint minimization problem of the 
instant method 10. 
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CLAIMS 

What is claimed is: 

1. In a method of digital image enhancement of a multidimensional digital image, said 
image being represented by a matrix [d] comprising image parameters, wherein said matrix 
[d] is mathematically manipulated to solve a linear ill-posed problem to reduce blurring, the 
improvement comprising: imposing a constraint on a reconstructed image matrix, said con- 
straint being based upon minimization of the area where strong variations and discontinuities 
between said image parameters occur. 

2. In the method of claim 1, the improvement further comprising: implementing said 
constraint in the form of weights imposed upon said image parameters. 

3. In the method of claim 2, the improvement further comprising: imposing penalization 
upon said image parameters, thereby to keep said said paramters within reasonable limits 
of variation. 

4. In the method of claim 3, the improvement further comprising: solving said ill-posed 
problem by means of an iterative loop using a programmed computer. 

5. In the method of claim 4, stopping said iterative loop when the norm of a differ- 
ence between the observed degraded image and a numerically predicted degraded image 
corresponding to an iteratively sharpened image reaches a tolerance value. 

6. A method of digital image enhancement of a multidimensional digital image, said 
image being represented by a matrix [d] comprising image parameters, comprising the steps 
of: 

a) initially restoring a digital image [mj by applying a transposed complex conjugated 
blurring operator, and an inverse gradient operator to the initial degraded digital image [d]; 

b) computing an inverse sharpening filter, thereby to minimize the area where strong 
image parameter variations and discontinuities occur; 

c) constructing a partially sharpened weighted image by applying said inverse sharpening 
filter; 

d) constructing an inverse filtered image by inverse filtering said partially sharpened 
weighted image using said inverse sharpening filter and said inverse gradient operator; 

e) checking the norm of a difference between the observed degraded image and a nu- 
merically predicted degraded image corresponding to said sharpened image; if said norm is 
equal to or less than a user defined tolerance value, then calculating the nonblurred image; 
otherwise, continuing to step f); 

f) undoing the results of loop steps comprising steps b), c), d), and e); and returning to 
step b). 

7. The method of claim 6 further including, subsequent to step d, and prior to step e), 
the additional step of imposing penalization to said inverse filtered image, thereby forcing 
the image parameters to be distributed within an interval bounded by a first upper value 
and a first lower value. 
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8. The method of claim 6 wherein said digital image is represented by a 2- dimensional 
matrix [d]. 

9. The method of claim 6 wherein said digital image is represented by a 3- dimensional 
matrix [d]. 

10. The method of claim 6 wherein one or more calculations comprised by steps of said 
method are performed using a programmed computer. 

11. The method of claim 6 wherein said inverse sharpening filter of step b is constructed 
using constraints comprising weights imposed upon image parameters. 

12. The method of claim 11 wherein said inverse sharpening filter comprises a minimum 
gradient support constraint. 

13. The method of claim 6 wherein step c further includes using conjugate gradient 
re- weighted optimization. 

14. The method according to claim 6, wherein: said step b)further comprises the steps 
of determining the spatial gradient of the partially sharpened image parameters and subse- 
quently deriving a variable weighting function accomplishing a minimum gradient support 
constraint based upon minimizing the area of a nonzero spatial gradient of the image pa- 
rameters. 

15. The method according to claim 14, wherein: said minimum gradient support con- 
straint is of the form ^^^^ = minimum. 

16. The method according to claim 6, wherein: said step c) futher includes the step of 
solving the following numerical formula for (n + l)-th iterations: 

= m£„ + <5m£ n = m£ n - fcj(m£ n ), 

wherein weighted conjugate gradient directions l(m% >n ) are determined by the expressions: 

and W m is the matrix of said sharpening filter, wherein the length of the step k n is determined 
by the line search to the minimum of the functional of the weighted image parameters 
P° (fri™ n - fc„J(m^ n )) . , and parameter (3 n is calculated by the formula: 




17. The method according to claim 6, wherein: said step d) further comprises recom- 
puting the real parameters of the ideal image from the weighted parameters m£ n+J at the 
(n + l)-th iteration by solving the following numerical formula: 

7fx n+1 = V-^^m^. 
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18. The method according to claim 7, wherein: said step imposing penalization comprises 
imposing the upper bound m" 6 (r) for the positive anomalous image parameter values, and 
the lower bound — m" (r) for the negative anomalous image parameter values, and discarding 
all the values outside these bounds during the entire image restoration process, thereby to 
keep the image parameters values m n (r) always distributed within the interval: 

m bg (r) - m l i (r) * m n (r) * m bg (r) + mf (r) . 
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Ideal seismic image of geological cross-section. 
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